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Chapter  1 


Introduction 


The  large  computational  requirement  for  the  simulation  of  a  wide  range  of  turbulence  [1],  acroaeous- 
tic  [2]  and  electromagnetic  [3]  phenomena  motivates  the  development  and  implementation  of  highly 
accurate  schemes.  In  direct  and  large-eddy  simulations  of  turbulence  for  example,  high-order  for 
this  purpose  fourth-order  accurate  or  higher  -  numerical  schemes  bring  some  previously  intractable 
problems  within  the  reach  of  modern  supercomputers.  An  overview  of  recent  efforts  may  be  derived 
from  a  study  of  Refs.  [4,  5,  6]. 

Although  the  complexity  of  problems  addressed  with  higher-order  schemes  has  increased  over  the 
past  several  years,  remarkably  few  studies  have  focused  on  wall-bounded  flows  around  geometrically 
complex  configurations.  Indeed,  higher-order  spatial  schemes  are  usually  coupled  with  explicit  time- 
integration  techniques.  This  is  usually  a  good  choice  in  situations  where  the  limiting  time-step  size 
is  dictated  by  physical  rather  than  numerical  constraints.  For  the  large  set  of  problems  constituted 
by  wall-bounded  flows  however,  the  stringent  mesh  resolution  requirements  near  walls  incur  a  severe 
numerical  time-step-size  limitation  which  can  be  alleviated  by  implicit  methods. 

The  objective  of  this  work  is  to  describe  the  incorporation  of  a  high-order  accurate  spatial  dif¬ 
ferencing  scheme  into  an  existing  implicit  flow  solver.  The  platform  chosen  for  implementation  is 
the  FDL3DI  code  which  is  the  primary  research  tool  of  the  CFD  group  at  Wright- Patterson  AFB 
(AFRL/VAAC).  The  basic  algorithm  employs  the  Beam- Warming  approximate  factorization  tech¬ 
nique  [7].  Several  enhancements  have  been  added  as  options  including  a  subiteration  technique 
and  an  efficient  diagonalized  procedure  [8].  Since  the  code  is  formulated  with  finite-differences, 
all  quantities  are  assumed  to  be  pointwise  in  nature.  Consequently,  the  formal  difficulties  encoun¬ 
tered  in  extending  finite- volume  approaches  to  higher-order  are  avoided.  However,  flux  and  metric 
conservation  issues  are  of  some  concern  and  are  addressed  elsewhere  [9]. 

Several  choices  arise  in  the  development  of  higher-order  schemes.  At  the  most  basic  level,  the 
formula  can  be  either  centered  or  upwinded.  While  each  has  its  advantages,  we  choose  centered 
formulas  because  of  their  nondiffusive  semidiscrete  error.  This  is  a  particularly  appropriate  choice 
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for  our  present  interests  which  encompass  relatively  low-speed  (i.c.,  subsonic)  and  thus  shock-free 
flows.  Extensions  to  include  shock-capturing  techniques  are  presently  under  development  and  will 
be  described  elsewhere.  For  a  fixed  order  of  accuracy,  centered  schemes  have  smaller  stencil  than 
upwind  schemes.  Even  within  centered  schemes,  additional  advantage  is  obtained  through  the  use 
of  "compact”  (or  Fade  type)  formulas  which  require  that  the  derivatives  be  computed  in  a  coupled 
fashion  along  an  entire  line  [4,  10].  With  this  approach,  greater  accuracy  is  obtained  with  fewer 
boundary  schemes  [10]. 

Compact  schemes  do  however  incur  a  moderate  increase  in  computational  expense  over  their  non¬ 
compact  counterparts.  In  order  to  limit  this  extra  effort.,  in  this  entire  work  formulas  are  restricted 
to  tridiagonal  systems.  With  this  .simplification,  a  particular  stencil  size  yields  two  orders  higher 
accuracy  than  an  explicit,  equivalent. 

The  formulas  required  to  treat,  various  aspects  of  the  solution  of  the  Navier-Stokes  equations  arc 
presented  in  Chapter  2.  In  developing  these  schemes,  the  approach  adopts  the  techniques  discussed 
by  Lele  [4].  The  principal  mathematical  tools  required  are  Fourier  analysis  and  Taylor  series  ap¬ 
proximations.  For  the  inviscid  terms,  formulas  are  required  to  evaluate  first  derivatives  of  the  metric 
quantities  and  the  fluxes.  The  maximum  stencil  size  for  this  clement  of  the  algorithm  is  chosen  to 
be  five  points,  the  highest,  scheme  so  obtained  being  of  sixth-order  accuracy.  The  coefficients  of  the 
five  point  scheme  can  be  adjusted  to  yield  lower  order  schemes  as  well.  Section  2.1  lists  the  various 
formulas  employed  in  this  work,  some  of  which  have  been  derived  previously  elsewhere  [10]  but  are 
reproduced  and  classified  within  the  framework  of  the  optimized  schemes  developed  in  Ref.  [1 1].  T  he 
interior  scheme  cannot  be  applied  at,  points  near  the  boundary  where  the  five  point,  stencil  protrudes 
the  domain.  These  special  formulas  have  not  been  standardized  to  the  same  extent  as  the  interior 
schemes  and  are  also  presented  in  Section  2.1. 

The  numerical  formation  of  the  viscous  terms  in  the  Navier-Stokes  equations  requires  calculation 
of  derivatives  of  the  components  of  the  shear  stress  tensor.  The  straightforward  approach  followed 
here  is  to  apply  the  formulas  of  Section  2.1  twice  in  succession.  However,  stability  considerations 
suggest  that  the  nonconservative  form  be  employed  together  with  formulas  which  compute  the  second 
derivatives  directly  (see  e.g.,  Ref.  [4]).  For  compact  schemes,  this  approach  is  expensive  either  in 
storage  or  in  number  of  operations  depending  upon  implementation.  An  alternate  strategy  is  to 
utilize  a  midpoint  interpolation  and  differentiate  sequence  to  evaluate  certain  viscous  terms.  These 
formulas  are  presented  in  Sections  2.2  and  2.3  respectively.  The  dual  formulas  to  compute  nodal 
derivatives  with  known  midpoint  values  are  described  in  Section  2.4. 

One  of  the  principal  problems  encountered  in  the  solution  of  the  Navier-Stokes  equations  with 
centered  schemes  is  the  appearance  of  numerical  instabilities,  typically  arising  near  boundaries  and 
in  regions  of  mesh  nonuniformities.  If  left  unchecked,  these  spurious  waves  contaminate  the  solution 
and  destroy  the  fidelity  of  the  solution.  A  common  method  to  suppress  such  instabilities  is  through 
artificial  dissipation  in  the  form  of  a  (small)  additive  damping  term  to  the  governing  equations  (e.g., 
Refs.  [12,  13]).  A  technique  of  similar  vintage  is  to  filter  the  solution  at  appropriate  intervals  in 
its  temporal  advancement  (e.g.,  Refs.  [14,  15]).  The  distinction  between  damping  and  filtering  is 
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relatively  subtle.  Ref.  [16]  notes  that  filtering  is  a  more  general  approach  not  restricted  to  hyperbolic 
equations.  In  Section  2.5,  tridiagonal  based  filters  of  up  to  tenth-order  are  presented.  In  each  formula, 
control  is  exercised  through  a  free  parameter  whose  range  and  impact  on  stability  are  investigated. 

The  implementation  of  the  various  formulas  into  the  FDL3DI  code  (to  yield  the  new  version 
FDL3DJCE)  is  described  in  Chapter  3.  An  additional  option  of  time-integration  has  also  been 
added  in  the  form  of  the  classical  fourth-order  Rungo-Kutta  method  (sec  c.g..,  R.cf.  [17]).  Finally, 
the  Appendix  presents  a  brief  description  of  the  input  parameters  to  the  FDL3DJCE  code  with 
particular  emphasis  on  those  which  are  either  new  or  whose  meaning  has  been  modified  from  the 
original  code.  Results  on  a  range  of  fluid  dynamic  problems  utilizing  simple  and  complex  mesh 
systems  can  be  found  in  Ref.  [9]. 
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Chapter  2 


Compact  Formulas 


Consider  a  1-D  mesh,  consisting  of  N  points  (or  nodes),  1,2,  . . 7  —  2,  i  -  1 ,  ?*  +  1 ,  7  +  2,  . . N  —  2, 

TV  —  1  (=  M),  N  as  shown  in  Fig.  2.1  (a).  Let  <f>  =  <£(*)  be  a  scalar  variable  whose  point-wise  values, 
<j)i  are  known  at  these  nodes.  We  assume  that  Xi+\  —  a:,-  =  J  i.c..  the  mesh  step-size  is  normalized 
to  unity.  For  body  conforming  meshes,  a  curvilinear  transformation  £  =  £(:r)  is  introduced  and  the 
same  formulas  are  then  employed  in  the  transformed  (£)  plane.  As  noted  earlier,  the  high-order 
method  requires  several  types  of  quantities,  formulas  for  which  are  now  obtained  in  succession:  a) 
first  derivatives  at  nodes,  b)  interpolation  of  quantities  from  nodes  to  midpoints,  c)  First  derivatives 
at  midpoints  in  terms  of  known  quantities  at  nodes,  d)  second  derivatives  at  nodes  from  known  first 
derivatives  at  midpoints  and  e)  filtering  formulas. 


2.1  First  Derivative  at  Nodes 

The  problem  is  to  utilize  the  known  <f>i  to  estimate  the  derivative,  <j>\  =  d<j>/dx\i  at  each  point  in  the 
mesh. 


2.1.1  Interior  Scheme 


At  interior  points,  a  centered  formula  is  employed: 


+  <t>\  +  CH<j)'i+1 


,<f>i  +  2  ~  <t>i-2  ,  <Pi+\  ~  4>i-l 

=  b - - - ha - - - 


(2.1) 


where  a,  a  and  b  are  constants  which  determine  the  spatial  properties  of  the  algorithm.  Note  that 
the  stencil  consists  of  five  points  as  shown  in  Fig.  2.1(b).  Up  to  sixth  order  accurate  schemes  can  be 
obtained  through  proper  choice  of  coefficients.  To  aid  in  this  procedure,  Taylor  series  approximations 
about  point  i  are  inserted  in  Eqn.  2.1  and  terms  of  various  orders  are  set  equal  to  zero.  This  gives 
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Figure  2.1:  (a)  Notation  for  1-D  Discretization,  (b)  Five-point  Stencil  for  First  Derivative  at  Interior 
Points,  (c)  Through  (f)  Stencils  for  First  Derivatives  at  Points  1,  2,  M  =  N  —  1  and  iV,  Respectively 
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Table  2.1:  Coefficients  for  Interior  Scheme.  OA=Ordcr  of  accuracy 


Scheme 

a 

a 

b 

Stencil 

OA 

Size 

ES 

0 

1 

0 

3 

2 

El, 

0 

4 

3 

-1 

3 

5 

4 

C4 

l 

A 

3 

2 

0 

3 

4 

GO 

1 

3 

\A 

9 

1 

0 

5 

6 

01 

0.351075 

1 .5673833 

0.1347667 

5 

4 

02 

0.381365 

1 .5875767 

0.1751533 

5 

4 

03 

0.347485 

1.5649900 

0.1299800 

5 

4 

01, 

0.370733 

1 .5804890 

0.1609770 

5 

4 

05 

0.430816 

1 .6205440 

0.2410880 

5 

4 

06 

0.376374 

1.5842493 

0.1684986 

5 

4 

07 

0.400218 

1.6001450 

0.2002910 

5 

4 

the  following  equations  (see  also  Ref.  [4]): 

0(/r)  :  1  -  a  +  2  q  -  6  =  0 

0{hA)  :  — a  -f  6  a  —  4  6  =  0  (2.2) 

0(h6):  -a+ 10  a- 16  6  =  0 

The  solution  of  these  equations  provides  values  of  a,  b  and  a.  Satisfaction  of  the  first  of  equation 
in  this  set  results  in  a  second  order  scheme.  If  the  second  equation  is  also  satisfied,  a  fourth  order 
scheme  results  and  the  unique  solution  of  all  three  equations  yields  a  sixth  order  scheme. 

Table  2.1  lists  several  standard  schemes  which  can  be  derived  by  appropriate  choice  of  coefficients. 
The  first  two,  E2  and  EJ  are  “explicit”  i.e.,  a  =  0  and  hence  the  derivatives  values  are  decoupled 
from  each  other.  For  the  remaining  schemes,  a  /  0  and  it  is  necessary  to  solve  a  tridiagonal  system. 
C4  is  the  original  fourth-order  compact  scheme  discussed  by  Hirsh  [10]  and  consists  of  a  three  point 
stencil.  C6 ,  described  also  in  Ref.  [4],  is  the  highest-order  scheme  obtainable  with  the  five-point 
formula  of  Eqn.  2.1.  A  semidiscrete  accuracy  analysis  in  the  context  of  the  wave  equation  has  been 
performed  in  Refs.  [4,  11].  Because  of  the  centered  stencil,  the  error  is  exclusively  dispersive.  The 
wave  propagation  speed  and  isotropy  characteristics  are  compared  with  the  exact  value  in  Figs.  2.2 
and  2.3  respectively  for  these  standard  schemes.  In  these  figures,  w  is  the  normalized  wave  number, 
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w  =  2irkh/ L  where  k  is  the  physical  wave  number  on  a  domain  of  length  L  and  h  is  the  grid  spacing. 
The  dramatic  improvement  in  the  compact  schemes  is  clear:  note  the  higher  accuracy  of  the  CJf 
scheme  compared  to  the  sixth-order  explicit  scheme  E6. 

Table  2.1  also  lists  a  values  for  several  fourth-order  optimized  schemes.  These  schemes  are  designed 
to  minimize  selected  error  quantities  over  various  wave  number  ranges  as  detailed  in  R.ef.  [18]. 
01  and  03  minimize  the  semidiscrete  isotropy  error  for  a  wave  spectrum  where  the  largest  wave  is 
resolved  with  4  and  |  intervals  respectively.  Equivalent,  schemes  which  minimize  dispersion  error  are 
designated  02  and  04 ,  respectively.  05  minimizes  the  dispersion  error  over  the  entire  range  of  wave 
numbers  up  to  2  points  per  wave.  However,  it  has  been  shown  in  Ref.  [1 1]  that  such  optimization  is 
counterproductive  because  the  absolute  error  is  prohibitively  large.  06  and  07  are  relevant  to  the 
fully-discrete  situation  where  the  time-integration  method  is  chosen  to  be  the  Runge-Kutta  classical 
fourth  order  method:  these  schemes  then  minimize  dispersion  error  at  GFL  numbers  /z  =  0.75  and 
1.0  respectively  for  a  wave  number  spectrum  resolved  with  four  or  more  points  for  every  component. 
Additional  discussion  on  these  optimized  schemes,  their  derivation  and  error  analyses  can  be  found 
in  Refs.  [11]. 

Special  formulas  are  required  at  points  1 ,  2,  —  1  and  N  where  the  stencil  of  Eqn.  2.1  protrudes 

the  domain.  These  formulas  constitute  numerical  (or  Phase  1)  conditions.  Physical  (or  Phase  II) 
conditions  arc  addressed  in  Section  2.6  and  in  Ref.  [9]. 

2.1.2  Boundary  Point  1 

In  order  to  maintain  the  tridiagonal  nature  of  the  scheme,  the  formula  employed  at  point  1  is: 

01  4*  —  ai0i  4-  b\<p2  4  c\(j>3  +  di<t>4  4  e\ 05  4  /i06  4  #i07  (2-3) 

The  stencil  is  shown  schematically  in  Fig.  2.1  (c).  Upon  inserting  Taylor  series  approximations 
about  point  1,  and  matching  coefficients  of  equal  order  terms,  a  sequence  of  equations  is  obtained 

whose  solution  yields  the  coefficients  which  are  listed  in  Table  2.2.  Again,  schemes  with  E  prefixes 

are  explicit  because  a i  =  0.  Note  that  the  scheme  C2  is  the  same  as  developed  in  Ref.  [19]. 

2.1.3  Boundary  Point  2 

The  general  formula  for  the  derivative  at  the  first  point  away  from  the  boundary  is: 

#2101  4  02  4  £*2203  —  <*201  4  ^202  4  C203  4  ^204  4  ^205  4  /206  4  #207  (2.4) 

Note  that  in  general,  both  sides  of  Eqn.  2.4  are  asymmetric  about  point  2.  Several  possibilities  arise 
from  different  utilization  of  </>[  in  Eqn.  2.4.  If  the  slope  at  point  1  is  treated  implicitly,  two  options 
arise: 

•  Option  A:  a^i  =  022  i1  0.  In  this  case,  the  left  hand  side  is  symmetric  about  point  2.  The 
coefficients  obtained  by  Taylor  series  coefficient  matching  are  presented  in  Table  2.3. 
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Figure  2.2:  Numerical  Versus  Theoretical  Wave  Number  Dispersion  Characteristics  of  Various 
Schemes.  normalized  wave  number,  PPW  =  points  per  wave 
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Figure  2.3:  Semidiscrete  Isotropy  Characteristics  of  Various  Schemes 
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Table  2.2:  Boundary  Coefficients  for  Point  1 


Scheme 

a  i 

6i 

Cl 

d, 

Cl 

/i 

9 1 

OA 

El 

0 

-1 

1 

0 

0 

0 

0 

0 

1 

E2 

0 

—  3 

2 

2 

-1 

2 

0 

0 

0 

0 

2 

E3 

0 

-11 

6 

3 

-3 

2 

1 

3 

0 

0 

0 

3 

E4 

0 

-25 

12 

4 

—3 

4 

3 

-1 

4 

0 

0 

4 

E5 

0 

-137 

60 

r> 

10 

3 

-5 

4 

1 

5 

0 

5 

E6 

0 

-40 

20 

6 

-15 

2 

20 

3 

-15 

4 

6 

5 

-1 

6 

6 

C2 

1 

-2 

2 

0 

0 

0 

0 

0 

2 

C3 

2 

-5 

2 

2 

1 

2 

0 

0 

0 

0 

3 

04 

3 

-17 

6 

3 

2 

3 

2 

- 1 

6 

0 

0 

0 

4 

05 

4 

-37 

12 

2 

3 

3 

—  2 
IT 

1 

12 

0 

0 

5 

C6 

5 

-197 

60 

—  5 

12 

5 

-5 

3 

5 

12 

-1 

20 

0 

6 

Table  2.3: 

Boundary  Coefficients  for  Point  2  with  Option  A: 

<*1 

=  a  2  ^ 

Scheme  <*21 

<*22 

G2 

62 

Co 

do 

e2 

/a 

92 

OA 

AC 4 

1 

4 

1 

4 

-3 

4 

0 

3 

4 

0 

0 

0 

0 

4 

A  C5 

3 

14 

3 

14 

-19 

28 

-5 

42 

6 

7 

-1 

14 

1 

84 

0 

0 

5 

A  C6 

2 

IT 

2 

IT 

-20 

33 

-35 

132 

34 

33 

—  7 

33 

2 

-1 

132 

0 

6 

10 


Tabic  2.4:  Boundary  Coefficients  for  Point  2  with  Option  B:  asi  ^  a->2  ^  0.  Note  BC4  scheme  is 
same  as  AC 4  because  this  formula  for  a  fourth-order  three-point  scheme  is  unique 


Scheme 

a  21 

O' 2  2 

a  2 

b  2 

Co 

do 

c2 

h 

92 

OA 

BC3 

5 

8 

-i 

8 

-3 

2 

3 

2 

0 

0 

0 

0 

0 

3 

BCJ, 

1 

4 

1 

4 

-3 

0 

3 

4 

0 

0 

0 

0 

4 

BC5 

1 

6 

1 

0 

-5  • 
9 

2 

1 

1 

18 

.0 

0 

0 

5 

BC6 

1 

8 

3 

4 

-43 

on 

-5 

6 

0 

8 

1 

6 

-1 

06 

0 

0 

6 

BC7 

1 

in 

1 

—  227 
600 

-13 

12 

6 

1 

3 

-1 

24 

1 

300 

0 

7 

BC8 

1 

5 

-70 

—  77 

55 

5 

-5 

1 

-1 

8 

12 

4 

240 

60 

48 

0 

48 

60 

720 

Table  2.5:  Boundary  Formulas  for  Point  2 

with  Option 

C: 

O' 2 1  — 

0,  022  7^ 

Scheme 

a  21 

C*22 

(In 

bo 

a  2 

do 

e2 

12 

92 

OA 

CC2 

0 

-1 

3 

—  2 
IT 

2 

3 

0 

0 

0 

0 

0 

2 

CCS 

0 

1 

2 

-1 

4 

-1 

5 

4 

0 

0 

0 

0 

3 

CC4 

0 

1 

_ 

-3 

2 

3 

2 

6 

0 

0 

0 

4 

CC5 

0 

3 

2 

-1 

8 

-11 

6 

3 

2 

1 

0 

-1 

24 

0 

0 

5 

CC6 

0 

2 

-1 

10 

-25 

12 

4 

3 

1 

-1 

6 

1 

60 

0 

6 

CC7 

0 

5 

-1 

-137 

25 

5 

-5 

1 

-1 

7 

2 

12 

60 

24 

3 

12 

12 

120 

#  Option  B:  <221  ^  022  j1  0.  Because  of  the  extra  degree  of  freedom,  for  the  same  stencil  as  in 
Option  A,  one  degree  higher  order  of  accuracy  is  obtained  as  shown  in  Table  2.4. 

In  most  algorithms,  the  solution  is  updated  only  at  interior  points  z.e.,  all  points  excluding  1  and  N 
where  the  values  <j>\  and  <j> jv  are  determined  from  the  physical  constraints  of  the  problem.  Thus,  the 
slopes  at  these  endpoints  are  not  required  and  it  is  possible  to  decouple  them  (but  not  the  pointwise 
values  <j)  1  and  <j>x)  from  the  rest  of  the  domain  by  setting  <*21  =  0.  Again,  two  possibilities  exist: 

•  Option  C:  c*2i  =  0,  <222  7^  0.  The  coefficients  of  this  implicit  formulation  are  listed  in  Table  2.5 


•  Option  D:  a2i  =  <*22  =  0.  In  this  case,  the  slope  at  point  2  is  computed  explicitly  with  the 
coefficients  listed  in  Table  2.6. 
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Table  2.6:  Boundary  Coefficients  for  Point  2  with  Option  D:  a 21  =  022  =  0 


Scheme 

Ooi 

O22 

a  2 

b 2 

(-2 

do 

h 

92 

OA 

DEI 

0 

0 

-1 

1 

0 

0 

0 

0 

0 

1 

DE2 

0 

0 

-1 

2 

0 

1 

2 

0 

0 

0 

0 

2 

DE3 

0 

0 

-1 

3 

2 

1 

-1 

6 

0 

0 

0 

3 

DEJ 

0 

0 

-1 

A 

—  5 

C> 

3 

2 

-1 

2 

1 

12 

0 

0 

4 

DE5 

0 

0 

-1 

5 

-13 

1  2 

2 

-1 

1 

3 

-1 

20 

0 

f> 

DE6 

0 

0 

- 1 

e 

—  77 
60 

n 

2 

-5 

3 

5 

e 

1 

A 

1 

30 

6 

2.1.4  Boundary  Point  M  —  N  —  1 

The  difference  scheme  at  the  point  away  from  the  right,  boundary  is  similar  to  that  derived  for  point 
2.  The  formula  chosen  is: 

QM 1  (t>  7yv  -  o  T  </>N  _  1  4  Ot  M  2&N  —  aM  0Ar  4  -  1  +  t'A/  - 2  4  d M  <f>  N  -  3  4  CM  <f>  N  -  A  4  fxi  <f>N- 5  4  -  f> 

(2.5) 

The  same  options  can  be  proposed  as  for  point  2.  Because  of  the  structural  relationship  between 
Eqs.  2.5  and  2.3,  Tables  2.3  through  2.6  arc  applicable  with  the  modifications  that  i)  a^i  =  »22< 
ii)  a a/ o  =  021,  and  iii)  the  signs  of  each  of  the  coefficients  a  through  g  are  reversed  i.e.,  clm  =  —0.2, 
bM  =  -62,  . . .. 

2.1.5  Boundary  Point  /V 

The  formula  for  the  boundary  point  TV  is: 

a/v^jv-i  4  =  a N<f>N  4-  +  cyv0;v_2  4  d/v</>jv_3  4  eN<f>N- a  4  /n^n-5  4  gN<}>N- 6  (2-6) 

Again,  because  of  the  structural  similarity  between  Eqns.  2.6  and  2.3,  Table  2.2  is  applicable  with 
i)  a/v  =  aq  and  ii)  the  signs  of  each  of  the  coefficients  a  through  g  are  reversed,  i.e.,  a/v  =  -fli, 
bN  -  — &! ,  . . .. 


2.2  Interpolation  Formulas  for  Midpoint  Values 

The  formation  of  some  viscous  terms  may  be  greatly  facilitated  by  the  use  of  function  and  derivative 
values  at  midpoints.  In  this  section,  the  function  values  are  obtained  at  midpoint  values  with  inter¬ 
polation  formulas  using  the  basic  procedure  of  Lele  [4].  The  notation  employed  and  the  schematic  of 
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Table  2.7:  Coefficients  for  Interior  Interpolation  Formula. 


Scheme 

a 

a 

b 

C 

OA 

C2 

Free 

J  4-  2o 

0 

0 

2 

C4 

Free 

9  ,  5a 

8  '  4 

-  l+6cr 

8 

0 

4 

C6 

Free 

5(1  5+1 4a) 

64 

—  25  + 126a 
128 

3—  10a 
128 

6 

C8 

5 

14 

25 

16 

5 

32 

-1 

224 

8 

the  stencil  are  sketched  in  Fig.  2.4.  We  focus  again  only  on  tridiagonal  schemes  with  up  to  8th-ordcr 
accuracy  i.e two  orders  higher  than  for  the  derivative  schemes  in  Section  2.1. 

2.2.1  Interior  Scheme 

The  basic  formula  for  interpolation  at  interior  points  is: 

*f  o^j+j  =  2  (&+i  +  +  2  (^+2  +  0*— l )  +  2  (0*+3  +  ^1-2)  (2.7) 

where  it  is  obvious  that  the  coefficients  0,  a  and  6  have  no  relation  to  those  of  Eqn.  2.1.  Matching 
of  Taylor  series  terms  gives  the  coefficients  listed  in  Table  2.7.  For  two  reasons,  a  is  retained  as  a 
free  parameter  for  all  but  the  highest-order  scheme  for  which  a  is  a  unique  nonzero  number.  The 
first,  advantage  concerns  notational  brevity:  explicit  schemes  can  be  easily  derived  by  setting  a  =  0 
while  the  equivalent  compact  scheme  can  be  obtained  by  setting  the  “outermost”  coefficient  to  zero 
-  thus  reducing  the  stencil  size  by  2  points  -  and  solving  for  a  unique  a.  For  example,  in  Table  2.7, 
an  explicit  6th  order  scheme,  E6,  can  be  derived  by  setting  a  —  0  in  the  row  labeled  C6 ,  thus  giving, 
a  =  75/64,  b  —  —25/128  and  c  =  3/128.  On  the  other  hand,  the  compact  6th  order  scheme  can  be 
obtained  by  setting  c  =  0  in  the  same  row  i.e.,  a  =  3/10  and  thus  a  =  3/2  and  b  =  1/10.  The  second 
reason  to  retain  a  as  a  free  parameter  is  that  it  may  then  be  potentially  employed  for  optimization 
purposes  similar  to  that  discussed  for  filtering  in  Section  2.5. 

2.2.2  Points  Near  Boundaries 

Point  | 

Formula 

$ |  +  —  a<j)  1  +  b<f)o  +  c<f) 3  4*  d<j) 4  4-  e<j)§  -1-  f<f)§  4-  g<j> 7  (2.8) 

Coefficients :  Table  2.8.  Note:  for  brevity,  the  subscripts  on  each  coefficient  have  been  omitted  in 
this  and  subsequent  formulas. 


13 


(a) 

Interior: 

(b) 


a 


a 


m 


i-2  i-1  i  i+1  i+2  i+3 

*  t  t  t  ♦  i 

-c  -b  -a  a  b  c 


N-2  M=  N 
N-1 


Implicit  stencil 


Explicit  stencil 


(c)  Point  3/2: 


1  a 


1  2  3  4  5  6  7 

i  i  i  k  i  k  . ♦ 

a  b  c  d  e  f  g 


(d)  Point  5/2: 


a 


a 


nm 


******* 


e  f 


Regular  point 


(e) 


Point  M-1/2: 


a  1  a 


ro 


N-6  N-5  N-4  N-3  N-2  N-1  N 

i  i  k  4— i 


g  f  e  d  c  b  a 


Regular  point 


(f) 


Point  N— 1/2: 


N-6  N-5  N-4  N-3  N-2  N-1  N 


t _ i _ i _ i _ i _ i _ f 

g  f  e  d  c  b  a 


Figure  2.4:  (a)  Notation  for  Midpoints  in  1-D  Discretization,  (b)  Stencil  for  First  Derivative  at 
Interior  Midpoints,  (c)  Through  (f)  Stencils  for  First  Derivatives  at  Points  |,  |,  M  =  N  —  |  and 
N  —  Respectively 
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Table  2.8:  Coefficients  for  Interpolation  Formulas  at  f  i.c.,  at  the  Midpoint  Between  Nodes  1  and  2 


Scheme 

a 

a 

b 

C 

d 

e 

/ 

g 

OA 

Cl 

Free 

1  4  o 

0 

0 

0 

0 

0 

0 

1 

C2 

Free 

1  _  a 

2  2 

1  +  2a 

2  '  2 

0 

0 

0 

0 

0 

2 

cs 

Free 

3  a 

8  8 

3  ,  3a 

4  4 

1  ,  3a 

8^8 

0 

0 

0 

0 

3 

Cl, 

Free 

5  a 

16  16 

15  ,  9a 

16  16 

5  i  9  a 

1 6  '  16 

1  a 

16  16 

0 

0 

0 

4 

C5 

Free 

35  5  cv 

128  128 

35  ,  15a 

32  '  32 

35  ,  45a 

64  "l_  64 

7  5a 

32  32 

5  i  3a 

1 28  '  128 

0 

0 

5 

CO 

Free 

63  7  a 

256  256 

315  ,  105a 
256  “r  256 

l5(-7-j-7a) 

128 

63- 35a 

128 

3(  —  1 5+ 7a) 
256 

7— 3a 

256 

0 

(i 

C7 

Free 

21(11  —or) 

63(1  l+3a) 

.  1  0  5  ( —  1 1  +9a) 

21  ( 1 1  —5a) 

9(  — 55+21  a) 

77—  27a 

7(  — 3+a) 

7 

1024 

512 

1 024 

256 

1024 

512 

1024 

C8 

1 1 

77 

693 

1 1 55 

—  77 

99 

-11 

7 

8 

3 

512 

256 

512 

128 

512 

256 

1 536 

Point  | 

Formula 

&<i>i  -f  d%  4  oid 7  =  a<^i  4  b<j> 2  4  c<j>3  4  d<f>4  4  e<t> 5  T  fds  4-  g<t> 7 


(2.9) 


Coefficients :  Table  2.9.  Note:  It  is  evident  from  Fig.  2.4  that  for  lower  even-order  schemes  (e.g., 
CV),  this  is  an  interior  point  and  the  formulas  are  the  same  as  in  Table  2.7. 

Point  M  -  \ 

Formula 


Old  M  —  %  40M-4  +  C^M  +  1  =  a<?>JV  4  ^N-l  4  cdN- 2  +  ddN~3  4  cdN-4  4  fdNS  (2.10) 
Coefficients:  Table  2.9.  Note:  the  coefficients  are  precisely  the  same  as  for  point  ~. 

Point  M  +  -  =  TV  —  - 
Formula 


OcdM-\  4  ^M+±  “  adN  4  bdN- 1  4  cdN-2  4  ddN-3  4  e<^>Ar-4  4  fdN-5  4  gdN-6  (2.11) 
Coefficients:  Table  2.8.  Note:  the  coefficients  are  precisely  the  same  as  for  point  §. 


2.3  Derivative  Formulas  at  Midpoints  from  Known  Nodal 
Values 

The  formation  of  the  stress  tensor  at  midpoints  also  requires  certain  derivative  values  at  midpoints. 
The  formulas  may  be  derived  in  essentially  the  same  manner  as  for  interpolation.  Since  the  node 
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Table  2.9:  Coefficients  for  Interpolation  Formula  at  Midpoint  % 


Scheme 

a 

a 

b 

C 

d 

C 

/ 

// 

OA 

Cl 

Free 

1  +-2n 

0 

0 

0 

0 

0 

0 

1 

C2 

Free 

=T  ~  » 

f  +  3o 

0 

0 

0 

0 

0 

2 

C3 

Free 

1  ,  3a 

8  '  4 

3  a 

4  2 

3  t  7a 

8  ’  4 

0 

0 

0 

0 

3 

C4 

Free 

1  ,  3a 

16  8 

9  t  5a 

16  ^  8 

9  ,  5a 

16'  8 

1  ,  3a 

16*  8 

0 

0 

0 

A 

05 

Free 

5  ,  19a 

128  '  6-1 

15  ,  15a 

32  *  16 

45  ,  5a 

64  32 

5  ,  11a 

32  '  16 

3  5a 

1 28  64 

0 

0 

5 

C6 

Free 

7  ,  33a 

256  “r  128 

105+290a 

1  5(7— 2a) 

-35+ 138a 

21 -70a 

-:i+l()a 

0 

a 

256 

128 

128 

256 

256 

\) 

C7 

Free 

7  (  — 3+3-la) 
1024 

7(27+94a) 

512 

31  5(3  — 2a) 
1024 

7(  —  1  5  + 58a) 
256 

l  89  —  670a 
1024 

—  27+98a 
512 

7  —  26a 
1024 

7 

08 

9 

21 

819 

945 

-21 

9 

-9 

1 

8 

38 

608 

1216 

1216 

608 

304 

1216 

1216 

'Iable  2.10:  Coefficients  for  Interior  Midpoint  Differentiation  Scheme.  OA=Order  of  accuracy 


Scheme 

a 

a 

b 

OA 

E2 

Free 

1 

0 

2 

CJ, 

Free 

9  3a 

8  4 

I  ,  n  a 

8  '  4 

4 

C6 

9 

62 

63 

62 

17 

62 

6 

derivatives  are  restricted  to  6th  order,  we  seek  only  up  to  sixth  order  accuracy  for  midpoint  deriva¬ 
tives  as  well. 


2.3.1  Interior  Scheme 

Formula 


+  <{>[+  i  +  +  l  =  a  (0i+ 1  -  <f>i)  +  g  (0i  +  2  “  0  i  —  1 )  (2.12) 

Coefficients  Table  2.10  Note:  Because  of  the  restriction  to  6th  order,  the  general  formula  does  not 
require  the  coefficient  c. 
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Table  2.11:  Coefficients  for  Midpoint  Differentiation  Scheme  at  Point  3/2.  OA^Order  of  accuracy. 


Scheme 

a 

a 

b 

C 

d 

C 

/ 

OA 

E2 

0 

-1 

1 

0 

0 

0 

0 

2 

C3 

23  ,  a 

24  24 

7— 9  a 

8 

1  ,  9a 

8  ‘  8 

—  1  —  a 

24 

3 

C4 

Free 

-22+a 

24 

1 7  9a 

24  8 

3-f  9a 

8 

-5- a 

24 

1 

24 

Co 

Free 

— 1689+71  a 

67-141  a 

143+ 207  a 

—  1  1 1  +a 

29  — 3  a 

—  71+9  a 

_ 

1920 

128 

192 

192 

128 

1920 

0 

C6 

62 

-10799 

-2713 

523 

-937 

25 

-3 

9 

17280 

384 

64 

1728 

384 

640 

0 

2.3.2  Points  Near  Boundaries 


Point  | 

Formula 


4-  =  a<j> i  +  h<j> 2  4-  4-  dd a  4-  ed>5  4-/06 


(2.13) 


Coefficients :  Table  2.11.  Note:  The  scheme  designated  (7/  is  not  really  compact  since  nr  cannot  be 
chosen  to  set  e  -  which  is  a  constant  -  to  zero. 


Points  |  and  M  —  \ 

For  all  orders  up  to  sixth,  these  are  interior  points.  Thus,  Eqn.  2.12  and  the  corresponding  coefficients 
of  Table  2.10  are  applicable  and  no  special  treatment  is  required.  Note  that  this  situation  differs 
from  the  case  for  interpolation  where  up  to  8th  order  formulas  were  derived. 

Point  M  4*  \ 

Formula 


0yv-i  4-  |  —  a<ps  4-  4-  c<j>N-2  4-  d<t>N-3  +  z4>n-a  4-  /0tv~ 5  (2.14) 

Coefficients:  May  be  derived  from  Table  2.11  by  reversing  the  signs  of  each  of  a  through  f  but  not 
of  a. 


2.4  Derivatives  Formulas  at  Nodes  from  Known  Midpoint 
Values 

The  midpoint  interpolated  and  derivative  values  obtained  from  Sections  2.2  and  2.3  can  be  employed 
to  form  composite  values,  in  this  case  certain  components  of  the  stress  tensor  as  outlined  later. 
It  is  necessary  to  then  differentiate  these  midpoint  composite  values  to  obtain  essentially  second 
derivatives  at  the  nodes. 
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Table  2.12:  Coefficients  for  Differentiation  at  Point  1  with  Known  Midpoint  Values 


Scheme 

Q 

a 

b 

c 

A 

e 

7 

OA 

C2 

Free 

-2-a 

3  +  o 

-1 

0 

0 

0 

2 

C3 

Free 

71  23  a 

24  24 

47.7a 

8  '  8 

—  3 1  +  a 

8 

23- a 

24 

0 

0 

3 

CJ, 

Free 

-93-22  a 

24 

229  i  17ft 

24  ^  24 

—  22+a 

24 

37  5  a 

8  24 

-22+a 

24 

0 

4 

C5 

Free 

—3043— 563a 
640 

5353+201  a 
384 

—3489+1 43  a 
192 

859 -37  a 

64 

—  204 1  +87  a 
384 

1  689  —  71  <•« 
1920 

5 

C6 

1627 

-1104667 

658913 

16343 

-6941 

15007 

-  10799 

(! 

62 

39680 

23808 

1 1  904 

3968 

23808 

119040 

u 

2.4.1  Interior  Scheme 

Formula 

° 1  +  4>i  +  $+i  =  r/  (<^i+i  “  ^i-i)  +  ^  -  0i-§)  (2.15) 

Coefficients :  Table  2.10.  Note:  At  interior  points,  this  is  the  dual  of  the  problem  of  Section  2.3  and 
thus  have  the  same  coefficients. 


2.4.2  Points  Near  Boundaries 

Point  1 

Formula 

<f>f]  4~  009  =  a03  -f-  b<f)$  4“  c*07  +  d<p9  +  e<p  jj.  4-  /0ii  (2.16) 

Coefficients :  Table  2.12 

Point  2 
Formula 

a<f>\  +  09  +  003  =  a0i  4-  60|  4-  c0z  4*  d<j> |  4-  e0jj.  4-  /0j^  (2.17) 

Coefficients:  Table  2.13. 

Point  M  —  N  —  1 
Formula 

Of0A/+ 1  +  0M  4-  Ot(j)'M_l  =  «0M  +  i  +  +  ^M-%  +  §  (2-18) 

Coefficients:  Obtained  from  Table  2.13  by  reversing  the  signs  of  each  of  a  through  f  but  not  of  o. 

Point  TV 

Formula 

<j)'N  +  o0yy _ j  =  a0yv_i  4-  &0/v-4  4-  c0N_§  4-  d(j) N_Z  +  e0^_|  4-  /0N-U  (2.19) 


18 


Table  2.13:  Coefficients  for  Differentiation  at.  Point  2  with  Known  Midpoint  Values 


Scheme 

a 

a 

b 

C 

d 

e 

/ 

OA 

E2 

0 

-1 

l 

0 

0 

0 

0 

2 

cs 

Free 

23  35  a 

24  12 

7  .  19a 

8  ‘  4 

1  11a 

8  4 

1  ,11a 

”  24  '  12 

0 

0 

3 

C4 

Free 

-11 -46  a 

12 

17  ,  101  a 

24  ^  12 

3  33  a 

8  4 

5  ( —  1 4-22  a ) 

24 

1  11a 

24  12 

0 

A 

C5 

Free 

—  1 689  — 9058a 

20 1+4930  a 

143- 3282 a 

-1  1  14-2578  a 

87- 2050  a 

-7i  +  1698« 

n. 

1920 

384 

192 

192 

384 

1920 

C6 

31 

5195 

4957 

119 

85 

119 

17 

6 

818 

'4908 

4908 

1227 

"  1227 

4908 

4908 

Coefficients:  Obtained  from  Table  2.12  by  reversing  the  signs  of  each  of  a  through  f  but  not  of  a  . 


2.5  Filter  Formulas 

The  interior  schemes  presented  in  earlier  sections  address  the  issue  of  accuracy.  An  equally  important 
consideration  is  that  of  stability.  For  high-order  schemes,  theoretical  analyses  of  stability  are  not 
straightforward,  the  principal  exceptions  being  the  simplest  cases  where  no  boundaries  are  present, 
the  governing  equations  are  linear  and  the  mesh  is  uniform  and  further  where  the  time-integration 
methods  are  explicit.  Practical  calculations  usually  do  not  satisfy  these  stringent  conditions  in 
several  ways.  Truncated  boundaries,  nonlinearities,  curvilinear  meshes  and,  in  the  case  of  wall- 
bounded  flows,  implicit  methods  of  time-integration  are  the  norm. 

As  noted  earlier,  special  schemes  are  necessary  at  points  near  the  boundaries  where  the  interior 
formulas  cannot  be  applied.  The  impact  of  the  formulas  employed  at  schemes  has  been  examined  in 
Refs.  [20,  21,  22]  and  the  citations  therein.  The  approach  is  typically  to  choose  a  well-posed  linear 
scalar  problem  discretized  on  a  uniform  discrete  mesh.  The  composite  scheme  is  then  subjected  to  an 
eigenvalue  analysis.  The  algebraic  complications  are  enormous  and  recourse  to  theories  connecting 
the  properties  of  semidiscrete  to  fully-discrete  schemes  are  often  invoked.  The  complications  are 
more  severe  for  systems  of  equations,  see  for  example  Ref.  [6]. 

Far  fewer  studies  address  the  imposition  of  the  second  or  physical  phase  of  the  boundary  condition 
implementation.  Typically,  these  include  physical  approximations.  For  example,  the  condition 
dp/dn  =  0,  where  n  refers  to  the  direction  normal  to  the  boundary,  is  imposed  at  solid  walls. 
The  condition  is  derived  from  boundary  layer  theory  and  although  it  is  used  extensively,  may  incur 
sizeable  error  near  points  of  separation  or  attachment.  Additionally,  the  actual  implementation  of 
this  condition  often  approximates  dp/dn  ~  dp/ dp  —  0  where  p  is  the  coordinate  line  emanating  from 
the  wall  and  may  not  point  along  the  normal.  The  impact  of  such  approximations  on  the  theoretical 
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stability  of  a  scheme  has  not  been  sufficiently  addressed  in  the  literature. 

The  quality  of  the  mesh  employed  also  plays  a  crucial  role  in  the  performance  of  the  scheme.  In 
Ref.  [23],  it  is  shown  that  central-difference  based  schemes  exhibit  spurious  reflections  at.  interfaces 
where  a  step-jump  in  the  spacing  is  encountered.  Indeed,  boundary  conditions  derived  for  uniform 
meshes  are  not  necessarily  stable  on  curvilinear  meshes:  Ref.  [24]  notes  one  such  instance. 

All  the  above  issues  introduce  uncertainty  int  o  t  he  applicability  of  tractable  theoretical  stability  anal¬ 
yses.  In  practice,  stability  is  usually  achieved  t  hrough  the  addit  ion  of  a  small  amount  of  damping  [13] 
or  less  commonly  through  filtering  which  forms  the  focus  of  this  study.  An  extensive  discussion  of 
the  basic  ideas  behind  filtering,  which  is  usually  applied  after  the  solution  vector  is  updated,  can 
be  found  in  R.cf.  [16]  and  the  citations  therein.  In  recent  years,  the  increased  use  of  very  high-order 
methods  has  encouraged  the  development  of  correspondingly  high-order  filters.  An  extensive  devel¬ 
opment.  of  explicit  filters,  i.c.  those  not  dependent  upon  the  solution  of  matrix  systems,  can  be  found 
in  Refs.  [16,  15]  while  the  methodology  to  derive  compact  schemes  has  been  outlined  by  Lele  [4].  In 
keeping  with  the  interior  algorithm,  this  work  utilizes  the  latter  approach. 

In  Ref.  [4],  a  series  of  optimized  tri-  and  pentadiagonal  compact  schemes  were  developed  for  up  to 
sixth-order  accuracy.  Several  were  designed  to  satisfy  specific  amplification  properties.  We  employ 
the  same  techniques  to  develop  a  different  set  of  filters  more  consistent  with  the  interior  scheme, 
restricting  attention  to  tridiagonal-based  schemes  of  up  to  10/7*  order.  Additionally,  the  variable 
aj  is  retained  as  a  free  parameter  in  order  to  provide  some  control  on  the  “degree”  of  filtering. 
With  these  choices,  a  2Nth  order  filter  has  a  stencil  of  2 N  +  1  points  and  the  stencil  is  wider  than 
that  in  Ref.  [4]  and  thus  require  derivation.  The  bounds  of  the  independent  parameter  aj  are 
established  and  the  performance  of  the  filter  is  characterized  for  various  values  within  these  bounds, 
in  effect  providing  guidelines  in  the  choice  of  this  parameter.  At  present,  the  conserved  quantities 
are  filtered  though  options  exist  to  apply  filtering  to  the  various  combinations  of  conserved  and 
primitive  variables. 


2.5.1  Interior  Scheme 


The  filtering  procedure  replaces  the  computed  (updated)  value  <f>  with  phi  obtained  from  the  following 
equation. 

Formula 

Otf(j)i-.  \  +  4>i  +  fl  =  £^-0““  (Ut+n  4*  Ui-n)  (2.20) 

Coefficients:  Table  2.14.  Note:  4>  are  the  filtered  values  of  </>.  The  spectral  function  (or  frequency 
response)  of  the  operator  is: 

£ff-0qncos (rvw) 


SF(w)  = 


(2.21) 


1  -h  2a/cos(tc) 

which  has  N  +  2  unknowns,  consisting  of  a/,  ao,  ai, . .  .  a;v.  To  obtain  the  coefficients,  we  first,  insist 
that  the  highest  frequency  mode  be  eliminated  by  enforcing  the  condition  SF{i r)  =  0  (see  Refs  [4, 
16]).  The  remaining  N  +  1  required  additional  equations  can  be  derived  by  matching  Taylor  series 
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Table  2.14:  Coefficients  for  Filter  Formula  at  Interior  Points,  aj  is  a  Free  Parameter  in  the  range 
0  <  a j  <  0.5. _ 


Scheme  a0  ai  a2  a3  a4  a5  OA 


F2 

i  +  aj 

h  +  aj 

0 

0 

0 

0 

2 

F4 

5  , 

8  '  A 

|  +  a/ 

+  £LL 

8  '  4 

0 

0 

0 

4 

F6 

11  , 

16  ■  8 

15  ,  17a/ 

32  '  16 

^3  ,  3a/ 

16  T  8 

1  Or/ 

32  16 

0 

0 

6 

F8 

93+70a  / 

128 

7+18a / 

16 

—  7+14a  / 

32 

1  or/ 

16  8 

1 28  '  64 

0 

8 

F10 

193+ 126a/ 

105+302a  / 

15(-l+2a/) 

■15(1  -  2a/) 

5(-l+2a/) 

1  —2a/ 

10 

256 

256 

6-1 

512 

256 

5 1  2 

coefficients  of  the  left  and  right  sides  expanded  about,  i .  Because  of  the  symmetric  nature  of  Eqn.  2.20, 
ST  is  real  and  the  filter  only  modifies  the  amplitude  of  each  wave  component.  Further,  only 
even  order  terms  persist  and  each  coefficient  match  yields  two  orders  of  accuracy.  The  property 
|5:F(u;)|  <  1  is  verified  a  posteriori.  For  concreteness,  we  choose  a  maximum  stencil  of  11  points 
with  seven  unknowns  (N  =  5).  While  the  seven  equations  arising  from  the  above  matching  procedure 
can  be  solved  uniquely  to  obtain  a  12th  order  filter,  the  formula  degenerates  to  an  identity,  aj  = 
1/2,  do  =  1,  ai  =  1/2,  an  =  0  for  n  >  1  and  ST( n)  =  §  is  indeterminate.  One  solution  is  to  suppress 
the  highest  coefficient  match  (thus  restricting  to  10th  order)  and  to  employ  instead  a  supplemental 
equation  by  assigning  ST(w)  at  some  specific  wave  number  as  in  Refs.  [4,  11].  Another  approach, 
as  in  Table  2.14,  is  to  retain  otj  as  a  free  parameter  which  provides  the  flexibility  of  an  explicit 
filter  subset  (obtained  by  setting  aj  to  zero).  Because  of  the  form  of  the  denominator  of  Eqn.  2.21, 
| or/ 1  <  0.5  .  The  spectral  properties  of  these  filters  have  been  examined  extensively  in  Ref.  [25]  and 
are  summarized  in  Figs.  2.5  and  2.6,  respectively.  Briefly,  the  damping  effect  decreases  with  increase 
in  order  of  accuracy  and  aj. 

2.5.2  Points  Near  Boundary 

For  an  11  point  stencil,  special  boundary  formulas  are  required  at  points  1, . . . ,  5  and  correspondingly 
at  TV  —  4, . . . ,  N  where  the  centered  filter  stencil  protrudes  the  boundary.  Two  approaches  to  treat 
these  points  are: 

•  At  any  boundary  point  i,  a  centered  scheme  of  order  2 i  —  2  can  be  easily  applied  with  the 
same  formulas  as  in  the  interior  (and  again  setting  all  free  parameters  with  the  exception  of  aj 
to  0).  Thus,  approaching  the  boundary  from  the  interior  this  strategy  systematically  reduces 
the  order  of  the  filter.  In  situations  where  the  mesh  is  refined  near  the  boundary,  there  is  no 
significant  loss  of  accuracy.  Further,  otj  can  be  considered  as  an  optimization  parameter  and 
at  values  close  to  0.5  can  provide  an  accurate  boundary  scheme  (see  Ref.  [9]  for  details). 
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Exact 


Figure  2.5:  Filtering  Effect  Variation  with  qj  for  8<h-order  Filter 
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Figure  2.6:  Filtering  Effect  Variation  with  Order  of  Accuracy  at  a/  =  0.25 
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•  One-sided  filter  formulas  can  be  incorporated  within  the  tridiagonal  structure  of  the  present 
scheme.  In  these,  for  general  values  of  ay  (still  satisfying  |ay|  <  0.5),  the  spectral  function  is 
complex  -  indicating  that  the  filter  may  introduce  artificial  dispersion  as  well  -  and  further,  the 
real  component  is  greater  than  unity  over  select,  wave  number  ranges  implying  that  certain  wave 
numbers  will  be  amplified.  Both,  the  imaginary  component  of  the  spectral  function  as  well  as 
the  amount  of  excess  over  unity  of  the  real  part  diminish  as  ay  approaches  0.5  and  thus,  a  value 
ofay  as  close  to  0.5  as  possible  from  practical  stability  observation  is  recommended.  However, 
it  is  also  important  to  note  that  since  these  schemes  arc  applied  at  only  a  few  points  in  the 
domain,  the  inaccuracies  introduced  at.  higher  wavenumbers  may  be  amply  compensated  by  the 
higher  accuracy  at  smaller  wavenumbers.  All  formulas  presented  below  guarantee  suppression 
of  the  odd-even  mode  as  for  the  int  erior  scheme.  Kach  formula  is  designated  F  Bx  y  where  x  is 
the  point,  where  the  formula  can  be  applied  and  y  is  the  order  of  accuracy. 

Point  5 

Formula 

a f<p4  4  05  -f  ay06  =  a<j>i  -f  602  4-  c^>3  4  d<f>A  +  c05  -f  /06  4  #07  4-  608  4  i<p9  4  j<j> m  4  (2.22) 

Coefficients  Table  2.15.  Note:  Only  two  schemes  are  listed  since  those  of  lower  order  can  be  obtained 

from  Table  2.14  by  considering  this  to  be  an  interior  point. 

Point  4 

Formula 

Qf<j> 3  4  <Pa  4  5  =  a<fi\  4  60 2  4  c03  4  d<j)A  4  e05  4  f<j>e  4  #0 7  4  h<b%  4  i<j> 9  4  j<f>  10  4  ^0i  1  (2.23) 

Coefficients:  Table  2.16. 

Point  3 

Formula 

Qtf(j) 2  4^3  4  tty0 4  =  a<j)  1  4  bfo  4  c03  4  d<j> 4  4  e0 5  4  f<t> 6  4  g<t> 7  4  h<j>%  4  i<j>$  4  j<j>  10  4  k<t>\  1  (2.24) 

Coefficients:  Table  2.17. 

Point  2 

Formula 

otf<j>  1  4  02  4  ay03  =  a<j>\  4  60 2  4  c03  4  d<j>A  4  e0s  4  f<f> 6  4  g<t> 7  4  h<p %  4  i<f> 9  4  j0io  4  k<f> n  (2.25) 

with  the  filter  coefficients  described  in  Table  2.18. 

Point  1 

Formula 

0i  4  07 02  =  fl0i  4  60 2  4  c03  4  d<f>A  4  e0 5  4  /0 6  4  g<f> 7  4  h(j>8  4  f09  4  j0io  4  &0n  (2.26) 
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Table  2.15:  Coefficients  for  Boundary  Filter  Formulas  at  Point  5.  aj  is  a  free  parameter 
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Table  2.18:  Coefficients  for  Boundary  Filter  Formulas  at  Point  2 
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Coefficients:  Table  2.19. 

The  filter  formulas  for  the  right  boundary  may  be  derived  in  similar  fashion  and  indeed,  the  same 
coefficients  are  obtained.  For  completeness,  these  are  listed  below. 

Point  N-4 

Formula 

<*/0/V-5  +  07V-4  4-  O/07V-3  =  G07V  4  607V- \  4  C0yv_ 2  4  d0yv- 3  4*  e<}>N-4  4  /07V- 5  4*  007V- 6  T 

/*07V-7  4  i<t>N-8  4-  /07V-9  4-  &0/V-1O  (2.27) 

Coefficients :  Table  2.15. 

Point  N-3 

Formula 

<yj<j>N-4  4-  0tv-3  4-  «/0tv- 2  =  «0tv  4-  60/v_i  4-  4-  d<j)N-3  4-  c0/v_4  4-  /0tv-5  4-  g<f>N-G  4 

h<i>N-7  4  *0tv-«  4-  j<t>N-9  +  A*0tv-io  (2.28) 

Coefficients :  Table  2.16. 

Point  N-2 

Formula 

a/0N- 3  4  0tv-2  4  a/^TV-i  =  <i0tv  4  60tv-i  4  c0tv_ 2  4  </0tv-3  4  e0/v_4  4  /07v- 5  4  g<t>N-e  4 
^07V~7  4  207V- 8  4  /07V-9  4  /u07V— 10  (2.29) 

Coefficients:  Table  2.17. 

Point  N-l 

Formula 

& f  07V  — 2  4  07V—  1  4  O/07V  =  a07V  4  ^07V—  1  4  C07V-2  4  d0TV_ 3  4  e0 7V-4  4  /07V  — 5  4  00 7V-6  4 
/?.07V-7  4  ?‘0/v-8  4  /07V  — 9  4  &07V-1O  (2.30) 

Coefficients :  Table  2.18. 

Point  N 

Formula 

O707V-1  4  07V  =  a0TV  4  60tv-i  4  c0tv_2  4  d<f>N_3  4  e07v-4  4  /0tv-5  4  00tv-6  4 

^0 TV— 7  4  207V  — 8  4  /07V  — 9  4  Ar07V—  i 0  (2.31) 

Coefficients:  Table  2.19. 

2.6  Boundary  Formulas  for  Neumann  Conditions 

Dirichlet  boundary  conditions  at  the  endpoints  1  and  TV  for  the  physical  (or  Phase  II)  conditions 
are  straightforward  to  implement  in  the  present  finite-difference  scheme.  Neumann  conditions,  such 
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Table  2.20:  Coefficients  for  Specification  of  06 /Ox  =  0 


OA 

a 

6 

a 

d 

e 

/ 

A 

6 

360 

-450 

400 

-225 

72 

-10 

147 

5 

300 

-300 

200 

-75 

12 

0 

137 

4 

48 

-36 

16 

-3 

0 

0 

25 

3 

18 

-9 

2 

0 

0 

0 

11 

2 

4 

-1 

0 

0 

0 

0 

3 

1 

1 

0 

0 

0 

0 

0 

1 

as  the  zero  pressure  gradient  condition,  require  higher  order  formulas.  The  following  formulas  sets 
the  value  of  <j>\  in  terms  of  interior  points  in  order  to  enforce  d<j>/dn  =  0. 


Point  1 

Formula 


Coefficients :  Table  2.20. 


4>\ 


a<f> 2  4-  b<j> 3  4  c<f>4  4  d6 5  4  c6e  4  67 

- 


Point  N 

Formula 


6n  = 


afiN- 1  4  b6 k-2  4  d6 n -3  4  e6N-4  4  /<ftyv-5  4  q6n- 6 

A 


Coefficients :  Table  2.20 


(2.32) 


(2.33) 


Chapter  3 


Implementation  for  Navier-Stokes 
Equations 


The  FDL3DI  code  was  chosen  to  implement  the  formulas  described  in  the  previous  chapters.  Several 
versions  of  this  code  exist,  some  including  capabilities  for  turbulent  simulations,  aeroelast.ic  analysis 
and  with  overset  grid  generality.  Some  of  these  versions  have  also  been  parallelized.  To  provide 
a  test  bed  to  investigate  the  properties  of  the  previously  developed  interior  differencing  and  filter 
schemes,  an  elementary  version  was  chosen  with  the  following  properties: 

•  Inviscid  fluxes:  Second-order  explicit  centered  or  Roe’s  flux-difference  split  scheme  with  all 
MUSCL-based  options. 

•  Implicit  time  integration  through  the  original  Beam- Warming  approximately  factored  algo¬ 
rithm  [26]  with  or  without  the  more  economical  diagonal  procedure  of  Pulliam  et  ai  [8]. 

•  Subiteration  capability  to  reduce  errors  introduced  by  linearization,  approximate  factorization 
and  explicit  boundary  condition  implementation. 

To  this  were  added  the  above  compact  differencing  and  filtering  capabilities  with  the  additional 
option  of  explicit  time-integration  with  the  classical  Runge-Kutta  scheme  in  the  low-storage  form 
described  in  Ref.  [27].  The  resulting  version  of  the  code  has  been  designated  FDL3DJCE. 


3.1  Governing  Equations 

The  code  solves  the  full,  unsteady,  three-dimensional  Navier-Stokes  equations  [28].  With  £,77,  f  and 
t  representing  a  general  curvilinear  transformation,  and  normalizing  with  the  quantities  by  poo,  600, 
a  representative  length  L,  and  the  free  stream  molecular  viscosity  p,  the  equations  may  be  written 
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in  vector  form  as: 


8  fU\  dF,  dGj  dfh  _  1  ,dFv  dGv  dH„ 
'dt\j)+dZ  +  dr1  +  dC  ~  Re[  d£  +  dr,  +  ^ 

where 


p 

pv 

Pv 

pW 

pu 

•  M 

puU  +iTp 

,  Gj  —  — 

puV  +  Tjrp 

•  «'  =  7 

pllW  +  Cr  P 

pv 

PVV  +iyP 

pvV  +  IjyP 

pvW  +  ( yP 

pw 

pwU  +i:P 

pwV  +  1];P 

pwW  +  C,zp 

.  pEt  . 

pE,U  +  pU 

pE,  V  +  pV 

pE,  W  +  pW 

0 

0 

0 

F  -I 
~  J 

tr.Rl 

1 

i 

Gji  i 

ixja 

Gja 

’  6’’' =  7 

1)xJi‘A 

,  it,,  =  - 

Gji-i 

Gja 

£r,b{ 

1)xJ>i 

Cr , 

U ,  V  and  W  are  cont.ravariant  components  of  velocity,  defined  by: 

U  -  it  +  £r«  +  iyV  +  i:W  =it  +0 

V  =  T),  +  T)rll  +  TjyV  +  T]:W  =  1],  +  V 

W  =  C t  +  Cr«  +  CyV  +  (zW  =  Cl  +  W 
E,  =  e  +  i(tr  +  v~  +  ur) 

Using  the  compact  notation  i  =  1 . .  .3  to  represent  the  x,  y  and  z  coordinates  respectively  and 
similarly  £,•  for  the  shear  stress  tensor,  r  can  be  written  as 

djkduj  ,  djk  duj ,  2  dji  duk 

Tij  ~  tl(dXj  dik  +  dxi  dik  ’  yij  dxk  dil 


while 


b,  —  UjT,j  T 


dii  dT 


(7  -  1  )PrMl  dxi  di, 

Note  that  the  pressure  is  normalized  by  the  quantity  PooU^-  The  equation  of  state  for  a  perfect 
gas  is  assumed,  the  Prandtl  number  is  fixed  at  0.72  and  the  molecular  viscosity  is  obtained  from 
Sutherland’s  law. 


3.2  Implementation 

3.2.1  Metrics  and  Inviscid  Fluxes 

The  metrics  are  evaluated  in  a  straightforward  fashion  by  utilizing  the  formulas  of  Chapter  2. 
Following  the  recommendations  of  Refs.  [29,  25],  it  is  recommended  that  the  same  equations  be 
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utilized  for  metrics  as  for  the  fluxes,  since  this  results  in  an  error  cancelation  effect.  The  derivatives 
of  the  in  viscid  fluxes  are  obtained  by  first,  forming  these  fluxes  at  the  nodes  and  subsequently 
employing  the  formulas  of  Section  2.1.  For  computational  efficiency,  the  derivatives  in  an  entire 
plane  are  obtained  in  a  vectorized  fashion.  To  maintain  the  overall  structure  of  the  original  code, 
the  and  values  are  obtained  in  a  tj  plane  by  vectorizing  in  the  (  and  £  directions  respectively. 
Similarly,  ^  are  formed  in  £  planes  by  vectorizing  in  £  . 


3.2.2  Viscous  Fluxes 


The  formation  of  the  viscous  fluxes  is  more  complicated  since  second  derivatives  are  involved.  In 
conservative  form,  and  with  x,  y  and  z  represented  by  X{ ,  i  =  1,2,3  respectively,  these  terms  have 
the  form  where  y  is  the  molecular  viscosity  and  ,s  denotes  one  of  the  primitive  variables 

u,  v ,  tv,  T .  At  present,  these  terms  are  formed  by  succesive  application  of  the  first  derivative  operators 
of  Section  2.1.  To  avoid  repeated  computations  of  derivatives,  the  various  terms  arc  formed  in 
different  subroutines  which  follow,  with  minor  exception,  the  original  layout  of  the  FDL3DJ  code. 
The  specific  routines  are  labeled  VCMPXZRHS  and  VCMPYRHS  with  some  terms  being  computed 
in  the  subroutine  VCROSX. 


It  has  been  suggested  that  the  successive  application  of  two  first  derivative  operators  can  lead  to 
instabilities  since  the  damping  mechanism  for  odd-even  modes  may  be  absent.  Of  course,  filtering 
guarantees  that  this  mode  is  suppressed  and  indeed  no  evidence  of  such  an  instability  was  found 
in  the  extensive  tests  reported  in  Ref.  [9].  The  suggested  alternative  to  the  successive  application 
of  first  difference  operators  is  to  write  the  chain-rule  form  of  the  viscous  terms  which  can  then  be 
directly  computed  with  formulas  for  second  derivatives.  We  note  however  that  this  is  an  expensive 
approach  either  in  CPU  or  memory  requirement  depending  upon  implementation,  particularly  in 
curvilinear  coordinates. 

An  approach  of  intermediate  complexity  can  be  suggested  with  the  use  of  the  midpoint  formulas 
presented  previously  (Sections  2.2  through  2.4).  The  cross-derivative  terms  can  be  computed  as 
before.  On  the  other  hand,  the  “straight”  derivative  terms  such  as  for  example, 


dxi 


(«|r)  =  (<“') 


/ 


can  be  discretized  by  first  forming  the  product  of  the  midpoint  values  of  y,  (with  the  formulas 
of  Section  2.2)  and  s'  (formulas  of  Section  2.3).  These  terms  of  the  stress  tensor  can  then  be 
differentiated  to  obtain  values  at  the  nodes  with  the  dual  formulas  of  Section  2.4.  The  result 
degenerates  to  the  original  second  order  method  with  the  appropriate  choice  of  coefficients  (scheme 
E2  of  Table  2.1). 
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3.3  Fourth-order  Runge-Kutta  Method 


The  classical  fourth-order  four-stage  Runge-Kutta  method  has  been  incorporated  into  FDL3DJ0E. 
With  R  denoting  the  residual,  the  governing  equation  is: 

dU  (d{Fj  +  Fv)  ,  d(Gj  +  Gv)  d(Ih  +  Hv)\ 

dt  ~  \  dt  ^  di]  +  dc  J 

The  classical  four  stage  method  may  be  written  as  [17]: 

kQ  =  AtR(Uo)  kl=AtR(U\) 

k2  =  A  tR{(/2)  ka  =  AtR{lh)  (3.2) 

(J 71  "i" 1  “  U o  +  “(i’o  *f  2k  \  +  2kn  +  h) 

6 

where  (Jo  =  Un ,  Ui  -  (Jo  T  k0/2,  (F  =  F i  +  Ar j  /2,  (/;t  =  U‘>  +  k2.  The  scheme  is  implemented  in  the 
low  storage  form  described  in  Ref.  [27].  The  residual  R  may  be  computed  with  any  of  the  schemes 
outlined  above. 
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Appendix:  Input  variables  for 
FDL3D_ICE 


The  specification  of  input  variables  is  illustrated  with  a  sample  input,  file.  Since  the  details  of 
this  file  are  subject  to  change  on  a  frequent  basis,  the  description  is  simply  to  provide  a  flavor  for 
the  required  data,  and  various  steps  necessary  to  formulate  a  computation.  The  input  file  contains 
comments  either  as  complete  lines  or  as  trailing  characters  which  describe  the  variable  name  being 
read.  In  the  first  case,  the  line  is  read  by  the  code  i.e.,  it  is  a  necessary  record  though  its  contents  are 
unimportant.  In  the  latter  case,  the  data  is  not  read  at  all  and  such  trailing  alphabetical  characters 
may  be  deleted  if  desired. 

A  few  routines,  such  as  the  one  applying  boundary  conditions,  are  case-specific  and  must  be  modified 
for  each  configuration.  These  routines  are  identified  below  in  the  appropriate  context.  Each  input 
noted  below  must  be  specified,  even  if  it  is  not  used  by  the  particular  scheme  selected.  For  example, 
if  one  of  the  variants  of  the  Roe  scheme  is  chosen,  the  compact,  difference  stencil  is  obviously  not 
utilized  but  must  necessarily  be  included  in  the  input  file  since  it  is  read  (though  not  utilized). 

In  the  following,  lines  in  the  input  file  are  entered  in  typewriter  font. 

4  TEST_CASE 

Variable  TEST-CASE  determines  the  case  being  computed.  This  term  is  only  employed  in  SUB¬ 
ROUTINE  BNDRY,  which  must  be  written  by  the  user  for  each  case. 

0.1  XM1 

100.0  RE 

1 . 002  TW 

0.38  SI 

0 . 0  ALFA 

XMl  is  the  Mach  number  and  RE  is  the  Reynolds  number.  TW  denotes  the  wall  temperature  Tw , 
normalized  by  the  free  stream  temperature,  and  is  employed  only  in  SUBROUTINE  BNDRY.  The 
convention  recommended  is  that  a  negative  TW  signifies  an  adiabatic  wall.  51  is  the  molecular 
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viscosity  coefficient  for  Sutherlands  law  (51  =  198.6/Too(Rankine)).  ALFA  is  the  angle  of  attack 
and,  like  TW,  is  only  employed  in  the  user  provided  SUBROUTINE  BNDRY. 

99999  1.  0  NDTAU  CFL  IBETA 

NDTAU  is  the  number  of  steps  between  update  of  tirncstep  size  (calls  to  SUBROUTINE  TMSTEP) 
and  is  meaningful  only  if  IBETA  ^  0.  CFL  is  the  Oourant  number.  IBETA  controls  whether  the 
time-step  size  is  controlled  with  the  CFL  number  or  if  constant  t  ime-step  size  values  DTVIS,  DTFIX 
(below)  are  employed.  For  further  details,  SUBROUTINE  TMSTEP  should  be  consulted, 

0.01  0.01  DTVIS  DTFIX 


DTVIS  is  the  time-step  size  if  IBETA=0  and  the  minimum  time  step  size  in  the  entire  domain  if 
local  stepping  is  chosen.  If  subitcrations  are  on,  DTFIX  becomes  the  outer  time  step  size  while 
DTVIS  is  the  step  of  the  “inner”  iterations. 

3  IDMPFIL 


This  parameter  controls  filtering/damping  execut  ion  and  consolidates  several  coefficients  of  the  pre¬ 
vious  version  of  the  code  including  ICDAMP  and  1SDAMP  as  follows: 


IDMPFIL 


1  (formerly  ICDAMP=1,  constant  coefficient  damping,  now  obsolete) 

2  (formerly  ISDAMP=1),  scalar  spectra!  damping  controlled  by  ES4,ES2, 

FES4I,  FES2I,  OMGAV  and  SRCONST  as  before 

3  Use  filtering  after  each  subiteration  (Note  implicit  damping  continues  to  be 
controlled  by  ES4,ES2,  FES4I,  FES2D,  OMGAV  and  SRCONST  as  before) 

4  Use  filtering  after  all  subiterations  are  completed  (Note  implicit  damping 
continues  to  be  controlled  by  ES4,ES2,  FES4I,  FES2D,  OMGAV  and  SR¬ 
CONST  as  before) 


0.01  0.0  ES4  ES2 

2.0  1.0  1.0  25.  FES4I  FES2I  OMGAV  SRCONST 


These  are  damping  coefficient  parameters.  These  values  are  always  relevant  for  implicit  approximate 
factorization  part  of  the  algorithm,  provided  the  time-integration  scheme  is  implicit  (controlled  by 
IRUNGE  below). 

0.  PHI 


PHI  controls  the  order  of  accuracy  of  the  implicit  time  integration  algorithm. 

1111  IV10N  IV20N  IV30N  IVM0N 

These  parameters  control  the  computation  of  viscous  terms  as  in  the  previous  version  of  the  code: 
IV10N(IV20N,IV30N)  =  1  switches  on  thin  layer  terms  in  I(J,K)  directions  respectively.  IVMON=l 
switches  on  the  cross  derivative  terms  as  well.  For  Euler  calculations,  IV10N=IV20N=IV30N— 0. 
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0 


1 


IDIAG  NCONV 


IDIAG=1  invokes  the  diagonalized  option  of  Ref.  [8]  while  a  zero  value  reverts  to  the  original 
Beam-Warming  scheme.  NCONV  controls  the  number  of  steps  between  runtime  output  (see  SUB¬ 
ROUTINE  WRTHIST  or  SUBROUTINE  MONITR). 

10  0  400  IRUNGE  ISUBON  NSUBMX  INMAX 

IRUNGE  controls  time  integration  (0  for  Beam-Warming  approximate  factorization),  1  for  RK4. 
Note  that  if  IRUNGE=1  (z.e.,  an  explicit  scheme),  many  other  inputs  which  control  the  implicit 
algorithm  are  moot  ( e.g IDIAG,  ISUBON,  NSUBMX).  ISUBON=l  switches  on  subiterations  whose 
number  is  controlled  by  NSUBMX.  The  parameter  is  ignored  if  IRUNGE=1.  INMAX  is  total  number 
of  iterations  for  the  present  run. 

15  15  15  1  1  ISCHME  JSCHME  KSCHME  IVISC  I METRO 

ISCIIEME,  JSCHEME  and  KSCHEMH  (previously  IROE,  JROE,  KROE)  control  the  choice  of 
scheme  as  follows: 

0  Original  second-order  central  scheme 

1  First  order  Roe 

2  Fully-upwind  second  order  Roe  (using  MUSCL  with  k  =  —  1) 

3  Second  order  Roe  scheme  with  Fromm  reconstruction  (k  =  0) 

4  Third  order  upwind-biased  Roe  scheme  (k  =  1/3) 

5  Second  order  central  subset  of  the  Roe/MUSCL  scheme  (k  =  1) 

15  Compact  schemes  as  determined  with  the  input  variables  below 

For  the  Roe/MUSCL  variants,  the  limiter  may  be  varied  by  modifying  MODULE  LIMTRS  in  the 
Fortran-90  version  of  the  code  or  in  SUBROUTINES  XIRECON,ETRECON  and  ZTRECON  of  the 
Fortran-77  version.  IVISC  controls  accuracy  of  evaluation  of  viscous  terms.  The  value  0  chooses  the 
second  order  central  formulation  of  the  original  code,  while  1  computes  the  derivatives  in  a  compact 
fashion.  Notei)  the  value  of  IVISC  is  ignored  for  the  /(</,  I<)  direction  if  IV10N(IV*20N,IV30N)=0. 
Also,  ii)  different  schemes  can  be  employed  in  different  directions,  including  different  variations  of  the 
compact  scheme.  However,  for  any  given  direction,  the  viscous  and  inviscid  fluxes  are  computed  with 
the  same  scheme.  Finally,  IMETRC  determines  the  difference  formula  for  the  metrics.  IMETRC=0 
reverts  to  the  original  second-order  scheme,  1  chooses  the  compact  scheme  determined  specially  for 
the  metrics  (below),  while  2  causes  the  code  to  compute  the  metrics  in  each  direction  with  the  same 
formula  as  for  the  inviscid  and  viscous  fluxes  in  that  direction. 

5.0E-2  5.0E-2  5.0E-2  XICUT0FF  ETACUTOFF  ZETACUT0FF 

000  IXIISO  IETAIS0  IZETAIS0 

000000000  ICUT1  ICUT2  ICUT3  JCUT1  JCUT2  JCUT3  KCUT1  KCUT2  KCUT3 

These  parameters  control  the  entropy  cutoff  for  the  Roe  scheme  for  each  of  the  three  directions.  De¬ 
tailed  formulas  can  be  obtained  from  SUBROUTINE  XIROEFLUX.  For  example,  in  the  I  direction, 
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XICUTOFF  controls  the  magnitude  of  the  cutoff,  IXIISO,  0  or  1,  controls  the  choice  of  isotropic 
or  anisotropic  formulas  while  ICIJTI,  ICUT2  and  ICUT3  are  on/off  switches  for  the  linear  (u)  and 
nonlinear  (u  +  c  and  u  -  c)  eigenvalues  respectively. 

!  THIS  LINE  BLANK 

FILE  NAMES:  INPUTD  (FOR  RECORD  ONLY ) , RESTART , RESTORE , SAVETMP 
’vcSd.dat’  ,-vc3d.in*  ’vcSd.out*  ’vcSd.tmp* 

These  are  names  of  files:  vc3d.dat  is  the  input  data  file  file.  The  code  must  be  executed  for  example 
with  the  command  “a.out.  <  vc3cl.dat”.  vc3d.in  is  the  restart  file  (must  exist  at  the  start  of  the  job 
with  a  flow  field  -  see  PROGRAM  FDL3DJOE),  v<*3d.out  is  the  restore  file  (created  at  the  end  of 
the  job).  vc3d.tinp  is  restore  file  dumped  every  50  iterations  (sec  PROGRAM  FDL3DJCE). 

!  THIS  LINE  BLANK 

MOVIE  INFO:  IMOVIE,  MODMOVIE  IMV(MIN/MAX/INT) , JMV(MIN/MAX/INT) ,KMV(MIN/MAX/INT) 

1  25  444-1-1-1  -1,-1,-1 

MOVIE  FILE:  GRID,  ROOT  NAME 
’movie . grid ’  ’ vc3dmv . J 

These  parameters  control  the  output  of  files  which  can  be  combined  to  form  a  movie.  The  details 
arc  highly  case-specific  and  may  be  tailored  by  appropriate  modification  of  SUBROUTINE  MOVIE. 

!  THIS  LINE  BLANK 

INFORMATION  ABOUT  COMPACT  DIFFERENCING  SCHEME 

****grid**** 

C4-AC4-C4-AC4-C4 

****x  DIRECTION**** 

C4-AC4-C4-AC4-C4 

♦***Y  DIRECTION**** 

C4-AC4-C4-AC4-C4 

****Z  DIRECTION**** 

C4-AC4-C4-AC4-C4 

These  inputs  select  the  compact  scheme  for  the  grid,  /,  J  and  K  direction  respectively.  The  character 
variable  consisting  of  five  fields  designating  the  scheme  to  be  employed  at  points  1,  2,  interior,  N 
and  N,  respectively.  For  example, 

C4-CC4-C6-AC4-C3 

requests  the  fourth  order  compact  scheme  at  point  1  (Table  2.2),  the  decoupled  fourth  order  com¬ 
pact  scheme  at  point  2  (Table  2.5),  compact  sixth-order  in  the  interior,  the  symmetric  compact 
fourth  order  scheme  at  point  N  —  1  (Section  2.1.4)  and  the  third-order  compact  scheme  at  point  N 
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(Section  2.1.5).  Note:  For  efficiency,  if  the  interior  scheme  is  explicit,  the  tridiagonal  system  is  not 
solved.  Thus,  it  is  not  possible  to  combine  explicit,  interior  schemes  with  implicit  boundary  schemes. 

!  THIS  LINE  BLANK 

PERIODIC  BC  INFO:  IPERDC , JPERDC , KPERDC 
1  1  1 

If  1,  IPERDC,  JPERDC  and  KPERDC  implement  periodic  boundary  conditions  in  the  j  and 
k  directions  respectively.  For  compact  schemes,  this  requires  the  solution  of  periodic  tridiagonal 
systems.  Note  that  a  5  pt  overlap  is  employed  such  that  point  1  and  2  corresponds  to  IEN  D  —  4 
and  l END  —  3,  respectively  while  points  IEN D  —  1  and  IEN I)  correspond  to  4  and  5  IPERDC, 
JPERDC  and  KPERDC  only  affect  the  compact  differencing  and  filtering  algorithms.  For  the  orig¬ 
inal  scheme  (ISCHME/JSCIIME/KSCIIME=0),  it.  is  necessary  to  explicitly  set  periodic  conditions 
in  SUBROUTINE  BNDR.Y. 

!  THIS  LINE  BLANK 

RK4  INFORMATION 

KSTAGES ,  COEFFICIENTS  as  in  CODE 

4  1.  6.  1.  3.  1.  3.  1.  6.  1.  2.  1.  2.  1.  1.  1.  2.  0.  1.  1.  1. 

These  inputs  are  the  coefficients  of  the  RK4  scheme  and  should  not  be  modified. 

The  next  inputs  control  the  application  of  the  filter.  Note  that  if  IDMPFIL=2,  the  filter  is  not 
invoked  in  any  direction. 


!  THIS  LINE  BLANK 

FILTER  INFORMATION 
I GF I LTER , JGFI LTER , KGFILTER 
1  1  1 

INOFIL , JNOFIL ,KN0FIL 
1  1  1 

At  the  highest  level,  filtering  can  be  switched  on  every  N  iterations  in  the  I(J,K)  directions  by  choos¬ 
ing  IGFILTER(JGFILTER,KGFILTER)=N.  If  IGFILTER=4  (for  example),  the  filter  is  executed 
every  4th  iteration  during  which  the  filter  is  applied  INOFIL  number  of  times  in  succession.  The 
significance  of  JNOFIL  and  KNOFIL  is  similar. 

FOR  I  DIRECTION:  INTERIOR  FILTER  ORDER,  ALPHA,  N00PTPARM,  OPTPARM(l) , . . . ,0PTPARM(N00PTPARM) 
10  0.4999  0 

These  inputs  control  the  filter  formulas  employed  in  the  I  direction  for  “interior”  points.  The  first 
parameter  is  the  order  of  accuracy  while  the  second  corresponds  to  ay.  The  10th  order  filter  (the 
highest  considered  here)  requires  a  stencil  of  11  points.  If  a  lower  order  filter  is  specified,  it  is  possible 
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to  provide  free  parameters  e.g.,  for  a  stencil  of  1 1  points,  an  8 th  order  scheme  permits  one  parameter 
to  be  free.  In  Table  2.14,  these  free  variables  (NOOPTPARM  in  number)  are  assumed  to  be  zero 
since  their  effect  on  the  filter  has  not  yet  been  investigated.  A  nonzero  NOOPTPARM  requires  spec¬ 
ification  of  the  appropriate  number  of  free  parameters  (OPTPARM).  For  details,  SUBROUTINES 
FILCOEFS  and  BCFILCOEFS  should  be  consulted. 

For  an  11  point  stencil,  interior  points  are  those  at  least  6  points  away  from  the  boundary.  It  is  thus 
necessary  to  specify  the  formulation  at  the  5  points  near  the  boundary.  As  noted  earlier,  at  each 
of  these  points  the  specification  can  proceed  in  one  of  two  ways.  The  first  approach  takes  recourse 
to  the  fact  that  points  near  the  boundary  can  be  considered  as  interior  points  for  lower  orders  of 
accuracy.  Thus,  at  a  point  A;,  a  centered  2N  order  formula  is  easily  constructed  from  Table  2.7. 

In  the  second  approach,  one  of  the  biased  formulas  presented  in  Section  2.5.2  may  be  utilized,  'lo 
allow  for  each,  the  code  inputs  are  specified  for  the  /-direction  as  follows. 

Point  1:  ORDER,  ALPHA\_F,  NOOPT,  OPTPARMS  (IF  NOOPT  >  0) 

0  0.0  0 
Point  2 
2  0.4999  0 
Point  3 
4  0.499  0 
Point  4 
6  0.49  0 
Point  5 
8  0.49  0 

Only  the  first  two  inputs  -  order  of  accuracy  and  aj  value  -  are  relevant  for  the  options  described 
here.  In  this  example,  the  order  of  accuracy  is  reduced  on  approaching  the  boundary  so  that  the 
formulas  are  always  centered.  To  compensate  for  the  loss  of  accuracy,  aj  is  successively  increased 
to  provide  a  relatively  “gentle”  filter  even  near  the  boundaries.  If  the  order  of  accuracy  is  higher 
than  achievable  with  a  centered  formula,  one-sided  equations  with  the  coefficients  of  Tables  2.15 
through  2.19  are  employed.  As  noted  earlier,  such  formulas  may  amplify  and  disperse  waves  and 
should  be  utilized  with  some  caution.  Thus,  at  present,  it  is  recommended  that  the  first  approach 
be  employed.  For  greater  flexibility  to  implement  optimized  filters,  NOOPT  can  be  set  to  a  nonzero 
value  together  with  additional  parameters  (OPTPARM (I),  1=1,  NOOPT).  Since  these  values  are 
not  yet  standardized,  this  option  is  not  recommended. 

To  permit  independent  specification  of  the  filter  in  the  different  directions,  separate  inputs  similar 
to  those  for  the  I  direction  illustrated  above,  are  employed  for  the  J  and  K  directions,  respectively: 

FOR  J  DIRECTION:  INTERIOR  FILTER  ORDER,  ALPHA,  SPECIAL  PARAMETERS  (number  and  values) 
10  0.4999  0 

Point  1  Order,  gamma,  noopt,  optparms 

0  0.  0 
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Point  2 
6  0.4999  0 
Point  3 
6  0.499  0 
Point  4 
6  0.49  0 
Point  5 
8  0.49  0 

In  this  case,  note  that  one-sided  formulas  will  be  invoked  at  points  2  and  3.  The  specification  of  the 
filter  for  the  K  direction  is  similar. 

FOR  K  DIRECTION:  INTERIOR  FILTER  ORDER,  ALPHA,  SPECIAL  PARAMETERS  (number  and  values) 
10  0.49  0 

Point  1  Order,  gamma,  noopt,  optparms 

0  0.  0 

Point  2 

2  0.4999  0 

Point  3 

4  0.499  0 

Point  4 

6  0.49  0 

Point  5 

8  0.49  0 

The  next  inputs  control  the  specification  of  H  cuts  in  the  domain.  This  is  a  useful  option  to  define 
grid-aligned  objects  in  the  field  and  to  blank  out  their  interiors. 

!  THIS  LINE  BLANK 

MULTIPLE  BLOCKS  N0BLKS , ( ( I JKBLK (IDMY , JDMY) , JDMY=1 ,6) , IDMY=1 ,N0BLKS) 

1  7  24  1  19  36  37 

NOBLKS  refers  to  the  number  of  such  cuts.  For  each  cut,  six  parameters  must  be  specified  denoting 
the  lower  and  upper  bounds  for  /,  J  and  K  indices  respectively.  These  parameters  are  employed  to 
automatically  blank  out  interior  points  for  the  implicit  scheme,  as  well  as  for  the  Phase  I  specification 
of  boundary  conditions.  Phase  II  specification  is  done  in  the  user-written  SUBROUTINE  BNDRY. 
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